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Abstract 

A new and very general technique for simulating solid-fluid suspensions has 
been described in a previous paper (Part I); the most important feature of the 
new method is that the computational cost scales linearly with the number of 
particles. In this paper (Part II), extensive numerical tests of the method are 
described; for creeping flows, both with and without Brownian motion, and at 
finite Reynolds numbers. Hydrodynamic interactions, transport coefficients, 
and the short-time dynamics of random dispersions of up to 1024 colloidal 
particles have been simulated. 



1 



I. INTRODUCTION 



In a previous paper (Ladd, 1993b), Part I of this series, the theoretical foundations of 
a numerical method for simulating solid-fluid suspensions was described. By combining 
Newtonian dynamics of the solid particles with a lattice-Boltzmann model of the fluid, a 
flexible and efficient algorithm was constructed, whose computational cost scales linearly 
with the number of particles. In this paper, the focus will be on the implementation of the 
algorithm for a variety of flow problems, and the numerical results. 

In section II, the numerical method described in Part I will be briefly recapped. Simu- 
lations of creeping-flow hydrodynamic interactions are described in section III; results are 
presented for periodic arrays, containing one or two spheres per unit cell, and for random 
dispersions of spheres. Section IV describes simulations at finite Reynolds number; flows 
past cylinders and spheres are compared with finite-difference and finite-element computa- 
tions at Reynolds numbers up to 200. Section V describes time-dependent flows, before the 
hydrodynamic interactions are fully developed, and section VI deals with fluctuations and 
Brownian motion. The work is summarized in section VII. 



In the lattice-Boltzmann approximation, the state of the fluid system is characterized 
by the discretized one-particle velocity distribution function n 8 (r,t), which describes the 
number of particles at a particular node of the lattice r, at a time t, with the discrete 
velocity c 4 -. The hydrodynamic fields, mass density p } momentum density j = pu } and 
momentum flux II, are moments of this velocity distribution: 



The specific model used in this work has 18 different velocities corresponding to the near- 
neighbor and second-neighbor directions of a simple cubic lattice; a complete list of the 
velocities is given in Table I. All quantities in this paper are expressed in "lattice units", for 
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which the distance between nearest-neighbor nodes and the time for the particles to travel 
from node to node are both unity. 

The time evolution of the distribution functions n 8 - is described by a discrete analogue of 
the Boltzmann equation (Frisch et al., 1987), 

n t (r + a, t + l)= m(r, t) + A,-(r, t), (2.2) 

where A; is the change in n 8 - due to instantaneous molecular collisions at the lattice nodes. 
The post-collision distribution n 8 - + A; is propagated to the neighboring node lying in the 
direction of the velocity vector c 4 -. A detailed description of the collision process is given 
in Part I. In essence, the collision operator is evaluated by linearizing about an appropriate 
equilibrium distribution, constructed so that 

/> = ZX 9 > j = ZX' c «'> w ' q = ZX' c '' c '' = p 1 + p uu - ( 2 - 3 ) 

i i i 

The equation of state is p = pc 2 s , where c s = yJl/2 is the speed of sound. 

Molecular collisions conserve mass and momentum but change the nonequilibrium (or 
viscous) part of the momentum flux Tl neq = II — IP 9 . In our simulations only the shear 
modes survive the collision process so that the post-collision momentum flux II' is given by 

IT = IT 9 + (1 + A)(TT - IT 9 ), (2.4) 

where the overbar indicates the traceless part of the tensor. The parameter A controls the 
rate of relaxation of the stress tensor; it is related to the shear viscosity of the fluid, 

r, = -p(2/\ + l)/6, (2.5) 

and lies in the range —2 < A < 0. In the low- Reynolds-number limit, fluid inertia is 
neglected; in this case Eq. 2.4 can be simplified to 

IT = pi + (1 + A)TT. (2.6) 

The post-collision distribution is given in terms of the mass density, momentum density, and 
updated momentum flux II', 
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rii + Ai(n) = a c 'p + a{i a c ia + a^W^c^c^ + a c 3 '(Tl' aa - 3pc 2 s ) : (2.7) 

the coefficients af are model dependent; values for the 18 velocity model used in this work 
are given in Eq. 2.11 of Part I. 

The solid particles are defined by a boundary surface, of any size or shape, which cuts 
some of the links between lattice nodes (see Fig. 2 of paper I). The fluid particles moving along 
these links interact with the solid surface at boundary nodes placed halfway along the links. 
At each boundary node there are two incoming distributions n 8 (r,t + ) and n 8 /(r + c 8 , 
corresponding to velocities c 4 - and c 4 -/ (c 4 -/ = — c 8 ) parallel to the link connecting r and r + c 4 -; 
the notation n 8 (r,t + ) = n 8 (r,t) + A 8 (r,t) is used to indicate the post-collision distribution 
(Eq. 2.7). The velocity of the boundary node u b is determined by the solid particle velocity 
U, angular velocity f2, and centre of mass R, 

u 6 = U + fix(r + §c,--R). (2.8) 

By exchanging population density between n 8 - and n 8 / we can modify the local momentum 
density to match the velocity of the solid particle surface at the boundary node, without 
affecting either the mass density or the stress, which depend only on the sum n 8 - + The 
precise form for the boundary-node collision operator is 

n t (r + Ci,t + 1) = n t f(r + c t ,t + ) + 2a c { pu b ■ c 8 , 

n 8 /(r, t + l) = ra;(r, t + ) - 2a{ pu b ■ c t ; (2.9) 

As a result of the boundary-node interactions (Eq. 2.9), forces are exerted on the solid 
particles, 

f(r + \ci,t + \) = 2[n t (r,t + ) - n t f(r + c t ,t + ) - 2a{ pu b ■ c t ]c t ; (2.10) 

thus momentum is exchanged locally between the fluid and the solid particle, but the com- 
bined momentum of solid and fluid is conserved. The forces and torques on the solid particle, 
obtained by summing f (r + |c 8 ) and (r + |c 8 ) xf (r + |c 8 ) over all the boundary nodes asso- 
ciated with a particular particle, are then used to update the particle velocity and angular 
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velocity, according to the laws of Newtonian mechanics. The mass and moment of inertia 
of the particle are preassigned, depending on the assumed distribution of mass within the 
particle. Since the forces at the boundary nodes are generated at the half-integer time steps, 
the force f at the intermediate integer time, 

f(r + \a,t) = (l/2)[f(r + \a,t- \) + f(r + \a,t + §)], (2.11) 

is used to update the particle velocities every two timesteps: 

U(t + 1) = U(t- 1) + 2M _1 F(*), n(t + 1) = n(t - 1) + 2I" 1 • T(t). (2.12) 

The particle mass M and moment of inertia I are preassigned parameters which control the 
rate at which particles respond to the fluid flow; usually M and I are on the order of several 
thousand (in lattice units). Since the velocities vary slowly on the time scale of a lattice- 
Boltzmann cycle, the precise form for the velocity update is usually not too important; 
however, it is important to use time-smoothed forces and torques, as described in Eq. 3.18 
of Part I. 

In this section we have briefly summarized the basic implementation of the simulation 
method. A detailed discussion of the method and its underlying theoretical justification is 
given in Part I. Methodology associated with specific applications will be described in the 
appropriate sections, along with the numerical results. 

III. LOW REYNOLDS NUMBER HYDRODYNAMICS 

Low- Reynolds number hydrodynamics is applicable to an intermediate range of particle 
sizes. For instance, in a suspension of particles sedimenting in water under the influence of 
gravity, the appropriate size range is between approximately l//m and 1mm. If the particles 
are smaller than about 1/mi, Brownian motion is important and the effect of fluctuations 
in the fluid must be considered (see section VI); for particles larger than about 1mm, fluid 
inertia begins to have an effect (see section IV). In this section three different sets of problems 
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are examined; periodic arrays of spheres, hydrodynamic interactions between two spheres, 
and the properties of random dispersions of spheres as measured by the short-time transport 
coefficients. The simulations were run under conditions such that the inertial contributions 
to the momentum flux were completely ignored (see Eq. 2.6). In all cases, results are 
compared with accurate numerical solutions of the Stokes equations, 



determined by multipole moment expansions of the force density on the particle surfaces 



In this section we focus on the translational and rotational friction coefficients of a 
simple cubic lattice of spheres. The force on a steadily moving sphere is computed under 
conditions of vanishing volumetric flow rate (Batchelor, 1972), which is enforced in two 
different ways. In the first case a quasi-periodic system is used, as illustrated in Fig. 1. 
Each solid particle is given the same translational velocity, perpendicular to the boundary 
walls. Since the system is Galilean invariant, the measured drag force depends only on the 
velocity difference between the solid particles and the boundary walls. The most convenient 
choice is to set the boundary wall velocities to zero, so that the volumetric flow rate is zero. 
Then the translational friction coefficient is just the drag force on the central sphere divided 
by its velocity. It has been verified that these simulations are unconditionally stable for 
any particle velocity or flow rate (set by the boundary wall velocity), and that the friction 
coefficient is independent of velocity. For more than three cells, the drag force on the central 
sphere is independent of the number of cells in the quasi-periodic array, implying that three 
cells is enough to simulate a fully periodic system. 

An alternate method is to hold the particles stationary (by considering them to be 
infinitely massive) and applying a uniform force density throughout the fluid, to simulate a 



V . u = 0, Vp = r]V 2 u, 




(Ladd, 1988; Ladd, 1990). 



A. Periodic arrays 
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pressure gradient across the system. For example, if a constant increment Aj x is applied to 
the x-component of momentum at each node, the macroscopic pressure gradient is given by 
V^p = Aj x . The fluid velocity at each node is measured after half the force has been applied; 
it was found that this symmetrical procedure led to an identical result to that obtained with 
quasi-periodic systems. The measured drag force on the particle Fd is smoothed over two 
successive time steps (see Eq. 2.11), 

F D = ^(l/2)[^(r 6 , t - \) + f x (r b , * + §)]; (3.2) 

here r b denotes the location of a boundary node. The balance of forces at steady state (for 
instance during sedimentation) implies that the total (gravitational) force on the sphere is 
balanced by the sum of the drag force and the buoyancy force Fb = — VsAj x (Vs is the 
volume of the sphere); thus 

F D + F B = (V- V s )Aj x , (3.3) 

or equivalently Fd = W x p. At steady state, the simulations satisfy this identity precisely. 
The volumetric flow rate JJy is measured from the time-smoothed flow field (see Eq. 3.20 of 
paper I), 

Uv = V- 1 ]T(l/4)Mr 6 , t - 1) + 2u x (r b , t) + u x (r b , t + 1)]. (3.4) 

Results from simulations at low solids volume fraction (less than 10%) are shown in Table 
II; numerical values from both quasi-periodic and pressure driven simulations were indistin- 
guishable. 

Since the particle surfaces are discrete, it is not possible to determine an accurate value of 
the particle radius a priori. An effective hydrodynamic radius a can be computed from the 
drag coefficient <^ T = Fd/Uv of a simple cubic lattice of spheres at low solids concentration, 
using the asymptotic expression for the drag coefficient (Hasimoto, 1959); 

6tt77 1 2.837 4.19 2 27.4 5 



The first term on the right hand side is the Stokes result for an isolated sphere; the next 
three terms are the periodic corrections to the Stokes drag for a unit cell of length L. In 
Table II three different estimates of the radius of the solid particle are given. The first 
is the input radius a which defines the boundary surface of the sphere; all lattice nodes 
inside this surface comprise the solid particle. Obviously, there is a range of values of a 
which lead to the same object, for instance y/2 < a < y/3. Thus a reasonable and well 
defined measure of the size of the object is the radius of the equal volume sphere ay. The 
effective radius, determined from Eq. 3.5, is also shown. All results in Table II are for 
A = — 1 (Eq. 2.4), corresponding to a kinematic viscosity v = 1/6 (Eq. 2.5). It can be seen 
that a consistent hydrodynamic radius emerges from these simulations, confirming that the 
simulated drag coefficients for dilute periodic suspensions are in agreement with Eq. 3.5. It 
should be emphasized that the same values for the effective radii of the objects are used in all 
subsequent simulations; the hydrodynamic radius a is not treated as a variable parameter, 
although even better agreement between theory and simulation could be achieved by very 
small variations in a. However, the effective radius is not completely independent of the 
fluid viscosity. Over a useful range of viscosities (10~ 3 < v < 1), a varies by about 1 lattice 
spacing (Table III). For large enough objects, this variation could be ignored, but in most 
practical instances, the effective radius must be calibrated for each fluid viscosity. 

In order to test the simulation method at higher solids concentrations, the translational 
and rotational friction coefficients of a simple cubic lattice have been calculated over the 
entire density range. The rotational friction coefficient was determined by measuring the 
steady-state torque for a fixed rotational velocity. The results are shown for various size 
spheres in Fig. 2. The overall agreement is quite good; spheres of radius 5-10 lattice 
spacings give quantitative results over the whole concentration range. It is also possible to 
"tune" the particle shape to a certain extent. For instance, the rotational friction coefficient 
for a sphere of radius a = 8.47, corresponding to an input radius a = 8.5 is about 10% 
too large at high solids concentrations . However, for a slightly smaller sphere a = 8.44 
(a = 8.4), the rotational friction is only about 5% too large. This difference arises because 

8 



of the discrete lattice, which, in certain directions, blunts the surface of the sphere into flat 
plates which have different areas for the two different radius spheres. This difference in area 
has a much greater effect on the rotational friction than the translational friction, because 
the faces are in relative motion. In random dispersions, these differences are less severe than 
for the lattice arrangements, and no attempt has been made in this work to optimize the 
particle shape. 

B. Hydrodynamic interactions 

A more stringent test of the simulations occurs if there is relative motion between the 
solid particles. In such cases there are singular lubrication forces when the particles are 
close to contact; forces along the line of centers diverge as s _1 and forces perpendicular to 
the line of centers diverge as Ins -1 , where s = (R12 — 2a) /a is the gap between the particles 
relative to the radius. The simulations comprise a periodic unit cell, 2L X L X L, with 
spheres located at (L ± L/2 } L/2 } L/2). The two spheres move with opposite velocities u 
and — u; thus there is no net momentum flux. For all velocities, the drag force was found to 
be a linear and superposable function of u. In Fig. 3 the simulation results are compared 
with integral-equation solutions (as in Ladd, 1988) for an identical geometry, including an 
explicit and exact calculation of the lubrication forces (Bossis and Brady, 1987). Again the 
overall agreement is good; the simulations are accurate for particle separations of I lattice 
spacing and beyond. However, the divergence of the lubrication forces near contact is not 
reproduced; instead the friction coefficients asymptote to a value more or less the same 
as that found at separations of I lattice spacing. These results could be improved upon 
by taking explicit account of the lubrication forces, as in Stokesian dynamics and related 
methods (Bossis and Brady, 1987; Ladd, 1990); however lattice-Boltzmann simulations of 
the short-range friction coefficients are much more accurate than typical multipole-moment 
or boundary integral methods and corrections for lubrication are not really necessary. Bulk 
properties of a suspension of such particles, which sample all distances, will be less sensitive 
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to the short-range interactions than the single-configuration results reported here. It will be 
seen that accurate simulations of the transport properties of dense suspensions can achieved 
with quite small particles. 

C. Transport coefficients 

In this section hydrodynamic transport coefficients of equilibrium distributions of spheres 
are reported. Experimentally these results correspond to short times, before any significant 
changes take place in particle configuration as a result of the externally imposed flow. We 
have computed the permeability of fixed random arrays of spheres (K), the collective mo- 
bility (//) or sedimentation velocity, the short-time self diffusion coefficient (D s ), and the 
high-frequency viscosity (?/oo)- I n Fig. 4, results from lattice-Boltzmann simulations are 
compared with independent simulations based on multipole moment expansions of the Os- 
een equation (Ladd, 1990); these latter results, when corrected for finite size effects, are in 
excellent agreement with experiment. Here we compare results for small periodic systems 
containing 16 spheres per unit cell in all cases; thus the plots in Fig. 4 differ somewhat from 
those reported earlier (Ladd, 1990), because of finite size effects. Typically, we average over 
10-100 different configurations of spheres for each calculation; the individual configurations 
were generated with a standard hard-sphere Monte Carlo program (Allen and Tildesley, 
1987). Thus we have a relatively rapid method of checking the accuracy of the simulations 
as a function of volume fraction <f> and particle size a. The results are clustered around 
three typical concentrations; dilute (</> ~ 0.05), semi-dilute (</> ~ 0.25), and concentrated 
(<f> ~ 0.45). 

The permeability of a fixed array of spheres relates the volume- averaged velocity (Uy) 
of the fluid (Eq. 3.4) to the pressure gradient 

U v = -— Vp; (3.6) 
V 

the velocity field is averaged over the whole volume of the system, including the interior of 
the spheres where the fluid is at rest. Thus the permeability of a random array of spheres is 
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calculated in a similar fashion to the translational friction coefficient of an isolated sphere; 
each sphere in the sample is held fixed, and a uniform force density = Aj is applied to ev- 
ery node at every time step. Numerical results show that the lattice-Boltzmann simulations 
predict the permeability of random dispersions of spheres rather well. Particles with radii 
as small as 2.5 lattice spacings are just as accurate as Stokesian dynamics (Phillips et al., 
1988), although larger spheres (a > 5) are needed for accurate results at high concentrations. 

The collective mobility (or sedimentation velocity) and self- diffusion coefficient also con- 
verge rapidly with particle radius, as shown in Fig. 4. It can be seen that the lubrication 
forces, which contribute in an average way to the self- diffusion coefficient, are quantitatively 
accounted for by the simulations. The methodology for determining the collective and self- 
mobilities is as follows. An external force F ext is applied to one particle (for D s = ksT fi s ) 
or is divided equally among all N particles (for //); this external force is balanced by a force 
density Aj = —F ext /V which is applied to all nodes in the fluid: thus 

F — l^si ^ — V- \0. I ) 

" ext " ext 

Many of the results for D s were computed by an alternate and comparable method in which 
a force F ext is applied to one sphere and a balancing force —F ext /(N — 1) is applied to all 
the other spheres. By measuring the average particle velocity of the tagged sphere < U\ >, 
we then obtain fi s from the relation 

Qttt] < Ui > _ (fi - fi s ) _ fi s N //_ 

F ext N-l ~N-1 N-l 

To simulate the motion of particles under the action of an external force, the fixed particle 
constraint must be removed and the particle velocities updated according to Eq. 2.12. In 
these calculations the mass and inertia were set to between 2 and 4 times that of an equal 
volume of fluid. The total force F is given by the sum of hydrodynamic force and the 
external force; at steady-state, these forces balance and the particle velocities are constant. 

The code is unstable for small values of either the mass or the inertia, because the 
boundary conditions (Eq. 2.9) reflect most of the incoming momentum. From a molecular 
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point of view this corresponds to the assumption that the solid particle is much more massive 
than the fluid molecules, which is not invariably so. It is straightforward to establish an 
approximate stability criterion. For a quiescent fluid (i.e. ignoring variations in /;), the 
momentum transfer per unit area per time step is approximately — 2/H1&. Integrating over 
the surface of the sphere we find, from Eq. 2.12, that the change in particle velocity AU is 
given by 

AU ~ --^U, (3.9) 

pa 

where f3 = p s / p is the ratio of solid particle mass density to fluid mass density. A similar 
expression holds for the angular velocity; the numerical factor in this case is 10. Thus we 
expect stable solutions if /3a > 10; it has been verified empirically that this is indeed correct. 

The high-frequency viscosity for fixed configurations of spheres in an external shear flow 
has also been computed. We have not, as yet, been able to devise a homogeneous shear 
algorithm, analogous to the those used for molecular liquids (Hoover et al., 1980; Ladd, 
1984). Thus we again use a quasi-periodic system with a geometry similar to that shown in 
Fig. 1, but with many spheres in each unit cell. Of course the configuration of spheres in 
each cell is the same. A shear flow is set up by moving one of the walls at a fixed velocity, 
u y} parallel to the wall. The stress and velocity profile in the central cell are measured, and 
from this the viscosity can be determined; once again results for quasi-periodic systems of 
3 unit cells and 5 unit cells are essentially identical. The stress in the suspension includes 
contributions from the fluid stress modes, together with particle fluid interactions; these are 
computed in a similar way to the torques, by summing r&f (r;,) over the particle surface. As 
before the simulations manage to pick up the effects of lubrication quite well, although at 
the highest volume fraction the smaller sphere (a ~ 2.5) is inadequate. A drawback of the 
present scheme for imposing an external shear flow is that, at higher solids concentrations, 
the velocity gradient is non-uniform in the outer two cells, although it is quite uniform in 
the central cell (Fig. 5). Thus the time it takes for the system to reach steady state can 
be quite long, as many as 10 5 -10 6 time steps, during which time a steady velocity profile 
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is slowly realized. It is to be hoped that a better choice of initial conditions will speed 
up the approach to steady state considerably. However, the slow approach of the initial 
configuration to steady state will be less of a hindrance in dynamical simulations, which can 
use long time averages from a single initial configuration instead of ensemble averages over 
many initial configurations. 

D. Shared nodes 

A technical difficulty arises in the simulation of many-particle suspensions which has 
not yet been discussed. As shown in Fig. 2 of paper I, boundary nodes sit at the midway 
points between lattice nodes. Thus some fraction of the boundary nodes can be shared 
between two different particles, with no intervening fluid between them. The number of 
"shared nodes" is generally small compared with the total number of boundary nodes; at 
fixed volume fraction the ratio is proportional to a~ 3 . At the highest concentrations (45%) 
the proportions of shared nodes are 50%, 5%, and 1% for spheres of radius 2.5, 4.5, and 
8.5 respectively. At lower concentrations, the proportion of shared nodes declines rapidly; 
thus for radius 2.5 spheres the proportion of shared nodes at solids concentrations of 45%, 
25%, and 5% is 50%, 8%, and less than 1% respectively. Although the specific details of 
the shared nodes are unimportant for sufficiently large particles (a > 10), the convergence 
is markedly improved by a careful treatment of these situations. 

To handle the shared nodes correctly, the effects of the local velocity of both particles 
must be taken into account. Of course, these velocities are similar (at steady state) so the 
boundary node velocity can be taken to be the average of that computed from the velocities 
of each particle (c.f. Eq. 2.8) 

u 6 = l - [U t - + n t x(r b - K t ) + U, + fi jX (r 6 - R,)] . (3.10) 

Using this velocity, the fluid populations are updated (Eq. 2.9) and the force is computed 
(Eq. 2.10). This force is then divided evenly between the two particles. The numerical 
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results lor small particles (a < 5) are quite sensitive to the exact prescription lor dividing 
up the force; variations on the even division ol force (for instance assigning the force first to 
one sphere and then, on the next time step, to the other) can produce significantly different 
transport coefficients for small spheres at high solids concentrations. 

IV. FINITE REYNOLDS NUMBER FLOWS 

The simulation technique can be carried over, without modification, to flows at non-zero 
Reynolds numbers. However, the simplified form for the equilibrium distribution (Eq. 2.39 
of paper I) cannot be used; the appropriate distribution is given in Eq. 2.5 of paper I. Two 
steady-state flows, up to Reynolds number R e = 200 have been studied; flows past a column 
of cylinders and flows past a simple-cubic array of spheres. Our results are compared with 
finite difference and finite element solutions of the Navier- Stokes equations. 

A. Column of cylinders 

The problem of fluid flow past a circular cylinder has been studied by many authors. 
Since our simulation method is designed for fairly dense suspensions, its use of a regular 
mesh of lattice nodes makes it rather inefficient for studying flows past isolated objects. A 
reasonable compromise is a column of closely spaced cylinders (Fornberg, 1991), which can be 
simulated by a long narrow channel with periodic boundary conditions. Our simulations are 
restricted to a single geometry, namely a column of cylinders separated by 10 cylinder radii. 
Test calculations with a small cylinder (a = 2.1) showed that drag coefficients within 1% of 
the infinite channel length result could be obtained for a channel length L = 16W = 160a. 
The flow was driven by setting the velocity of a single row of nodes to the desired input 
velocity Uq\ the cylinder was placed 411^ = 40a downstream, leaving a maximum wake length 
of 1211^ = 120a. Simulations were run with cylinders of input radii a = 1.5, a = 4.5, and 
a = 9.5 with a kinematic viscosity v = 0.01. The effective hydrodynamic radii of the 
cylinders were determined to be 2.1, 5.1, and 10.1 respectively. The calibration of the 

14 



cylinder size was performed in an analogous way to that for spheres, by measuring the drag 
force of a dilute periodic array of cylinders at R e = 0. To a first approximation (in solids 
concentration) the friction coefficient £ T is given by the two-dimensional equivalent of Eq. 3.5 



2 1= 1_ 

A ^ 



Jo(ka) 



(4.1; 



k 

A is the area of the unit cell, k = 2mr / L is the wavevector, and J is a cylindrical Bessel 
function. 

A contour plot of the stream lines for the larger (a = 10.1) cylinder is shown in Fig. 6 
at various Reynolds numbers (R e = 2Uoa/p); the maximum Mach number (at R e = 100) 
is 0.07. The flow at R e = 100 is eventually unstable; however, the instability is sufficiently 
delayed that a quasi-steady state could be reached before its onset. A detailed comparison 
of the simulation results with Fornberg's finite difference code is shown in Table IV. 1 It can 
be seen that the lattice-Boltzmann simulations reproduce all the measured characteristics 
of the flow quantitatively, including the minimum in the stream function at the center of 
the recirculation zone. 



B. Simple-cubic lattice of spheres 

These simulations are similar to the calculations of the creeping-flow translational friction 
coefficient described in section IIIA. The sample contains a single stationary sphere of radius 
L/2 in a cubic box of length L with periodic boundary conditions; the flow is driven by a 
uniform force density, as before. The only difference is that the full non-linear form for the 



1 Fornberg's published results (Fornberg, 1991) begin at R e = 100. Because his code utilizes the 
symmetry present in steady flows, he obtains time independent solutions at all Reynolds numbers. 
Our simulations exhibit vortex shedding at higher Reynolds numbers; we can only obtain steady 
solutions up to R e = 100. Dr. Ramesh Natarajan kindly ran some intermediate Reynolds numbers 
for these comparisons. 
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equilibrium distribution is used, with a kinematic viscosity v = 0.01; the Mach number is 
then small, typically of the order of 10~ 2 . The results are reported in Table V. The input 
radii (a = L/2 — 0.2) were determined so that the effective radius a was as close as possible 
to L/2; the choice of input radius is consistent with the result reported in Table III for 
v = 1/96. 

It can be seen that the mean flow velocity JJy (or Reynolds number) for a fixed drag force 
(Fd = — Wp) converges quite rapidly with increasing particle radius (or unit cell length L). 
Thus it is surprising that the simulation results do not agree at all with the results of Lahbabi 
and Chang, 1985, who studied the identical problem via a truncated mode expansion of the 
Navier-Stokes equations. A comparison is shown in Fig. 7; the reduced pressure gradient 
Vpa 3 / pv 2 is plotted vs Reynolds number R e = 2LVa/ez/, where e = 1 — 7r/6 is the porosity 2 . 
At Reynolds numbers greater than about 40 the simulations deviate rapidly from Lahbabi 
and Chang's results; they differ by about a factor of two at R e = 200, with the simulations 
predicting a much smaller pressure drop. To investigate this discrepancy further, a finite 
element model of the problem (Tompson, 1983) has been run; these results, shown as open 
symbols in Fig. 7, are in excellent agreement with the lattice-Boltzmann simulations 3 . The 
lattice-Boltzmann/finite element results together suggest that the calculation of Lahbabi 
and Chang is incorrect. Moreover, they indicate that finite Reynolds number corrections 
to the Stokes drag in a regular arrays of spheres are considerably smaller than in random 
arrays; only about 20% at R e = 100. 



2 Note that Lahbabi and Chang use a different definition of R e , based on the particle radius rather 
than the diameter; thus their Reynolds number is smaller than ours by a factor of 2. 

3 The finite element code LAMFLOW used for these comparisons is based on a penalty-function 
formulation of the time-independent Navier-Stokes equations (Hughes et al., 1979). The calcula- 
tions used 148, 8-node prismatic elements to represent one quadrant of the periodic unit cell. The 
code was recently resurrected for this problem by its author, Dr. Andy Tompson 
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V. TIME-DEPENDENT HYDRODYNAMIC INTERACTIONS 



The classical description of colloidal particle dynamics is the Einstein- Smoluchowski 
equation, in which the details of the short time dynamics are ignored. The hydrodynamic 
interactions are assumed to be fully developed; in other words there is a complete separation 
of time scales between the dynamics of the fluid and the (diffusive) motion of the solid parti- 
cles. In reality, hydrodynamic interactions develop by the diffusion of fluid vorticity, with a 
characteristic time scale tj = a 2 jv. Lattice-gas/lattice-Boltzmann simulations methods are 
motivated by the observation that, by solving time-dependent Navier-Stokes equations, the 
correct hydrodynamic forces and torques can be computed from purely local interactions; 
thus the computational cost scales linearly with the system size. However, for low Reynolds 
number flows, an extra time scale separation between the solid and fluid motions must now 
be introduced, to allow time for the steady-state hydrodynamic interactions to develop. In 
practice this is not a severe constraint, as can be seen from the following considerations. 
The time scale characterizing the particle motion is t s = a/U } where U is a typical particle 
velocity; the ratio tj/t s is just the Reynolds number R e . It is known that the creeping-flow 
limit is generally equivalent to conditions in which R e < 1 in dense suspensions and R e < 0.1 
in dilute suspensions (Happel and Brenner, 1986); for a particle radius a = 5 and typical 
kinematic viscosity v = 1/6, this implies particle velocities less than 0.01. Thus a typical 
dynamical simulation (Ladd, 1993a), where the particles must travel distances on the order 
of 100 radii, would require about 10 5 time steps. Based on a particle radius a = 4.5, a com- 
plete sedimentation run would require an estimated total of 10 8 node updates per particle; 
a modern workstation (IBM 580 or HP 735) can update more than 10 5 nodes per second, so 
that a complete run would take about 15 minutes per particle. Moreover, the code can be 
readily executed, in parallel, on many processors at once. Thus simulations of 10 5 or more 
spheres should be feasible on a massively-parallel supercomputer. Even with the present 
code, systems of 128 spheres (see section VI) are simulated routinely (mainly on SPARC I 
and SPARC II computers), and in a few instances systems of 1024 spheres (on an IBM 580). 
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Although the original motivation for developing a method that simulated time-dependent 
hydrodynamic interactions was purely computational, the emergence of diffusing- wave spec- 
troscopy (Weitz et al., 1989) has initiated experimental studies into the motion of colloidal 
particles at very short times. Underlying the experimentally observed behavior (Zhu et al., 
1992; Kao et al., 1993) are both time and space dependent hydrodynamic interactions; with 
this new simulation technique we have a unique capability to determine these interactions 
ah initio (Ladd, 1993c). Since these experiments track the Brownian motion of the colloidal 
particles, a brief discussion of the relevant simulations will be deferred to section VI; an 
account of this work has already been reported (Ladd, 1993c). Here we focus on studies 
of the motion of isolated spheres, to verify the accuracy of the simulations in tracking the 
temporal and spatial development of the hydrodynamic forces (see also Fig. 5 of paper I). 

The unsteady motion of an isolated sphere has been worked out in some detail; the 
basic results for low Reynolds number (ignoring fluid inertia) can be found, for example, 
in Landau and Lifshitz, 1959. Recently this work has been extended to small but finite 
Reynolds number (Lovalenti et al., 1993); however, at present our results are limited to 
the time dependent R e = case only. We begin by investigating the motion of a sphere 
undergoing small amplitude rotational oscillations 0(t) = O cos(o;t); the in phase and out 
of phase components of the drag torque were measured to determine the modulus of the 
torque |T| and the phase lag ip (see Fig. 8). Since there is fluid both inside and outside the 
sphere, contributions to the drag torque from both regions are summed together to give the 
solid lines in Fig. 8. These expressions, in terms of the reduced frequency u* = coa 2 / ' v and 
the zero frequency drag coefficient £ R = 8irr]a 3 , are (Landau and Lifshitz, 1959) 



It can be seen (Fig. 8) that the simulated phase lag is inaccurate at sufficiently high frequen- 
cies, when the period of oscillation (2tv/u) is less than about 10 simulation time steps. High 
frequency motion can be accurately simulated either by increasing the size of the spheres or 
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by decreasing the viscosity of the fluid. However, reduced frequencies (ua 2 jv) larger than 
10 are unimportant in most instancies. 

At low frequencies, the drag forces and torques from the interior fluid reduce to the 
additional inertial drag of the oscillating mass of fluid; for rotational oscillations the interior 
torque becomes Ti nt = 8tt pa 5 icoO, / '15 (Eq. 5.1). Since the added inertia from the interior 
fluid can be readily accounted for (as will be seen shortly), only the deviations in the interior 
drag force from the inertial drag are significant; these deviations are proportional to cu 2 , to 
leading order. Thus the effects of the interior fluid can be minimized by using a smaller 
viscosity in the exterior fluid than in the interior fluid as shown in Fig. 9. Here the interior 
fluid viscosity is kept fixed at v = 1/6 whereas the exterior fluid viscosity is varied between 
1/24 and 1/96. It can be seen that for the same size sphere (a = 4.5) agreement between 
simulation and theory improves rapidly with increasing viscosity ratio (see also results in Fig. 
8 for a = 4.53); moreover deviations of the drag torque from the asymptotic form (exterior 
torque plus interior inertia) become small. However, it will be seen that this refinement is 
unnecessary in time-domain simulations, since these high frequencies play a very small role. 

In Fig. 10 the decay of an initial translational velocity U(0) or rotational velocity 0(0) 
of a sphere is compared with theoretical predictions based on an inverse Laplace transform 
of the frequency-dependent equations of motion (Hauge and Martin-L6f, 1973) 



The friction coefficients that appear in Eqs. 5.2 are the exterior friction coefficients 






(5.3) 



the internal drag forces are approximated as purely inertial mass; iloMj = Att pa 3 iuj /3 for 
the translational friction and iloIj = 8tt pa 5 ico / '15 for the rotational friction. 
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The decay of the translational velocity of a sphere can be expressed analytically (for 
example Hinch, 1975), but the rotational velocity has not been computed over the full 
time domain. Thus we follow Hauge and Martin-L6f, 1973 and convert the inverse Laplace 
transform to a definite integral (by contour integration) and then evaluate this integral by 
numerical quadrature to get the solid lines in Fig. 10. The mass and inertia of the solid 
particle comprise the sum of the assigned mass M or inertia /, used to update the particle 
velocities (Eq. 2.12), and the mass Mf = Airpa 3 /3 or inertia If = 8irpa 5 /15 of the interior 
fluid. Although there are no free parameters in the comparison of simulation with theory, 
the agreement is essentially perfect over the whole time domain. This supports the earlier 
assertion that the interior fluid contributes an inertial force, due to its extra mass, and very 
little else. 

It is appropriate at this point to comment on the normalization of the theoretical results; 
they are not normalized to the t = velocities U{0) and 0(0) but to something slightly 
less. There are two reasons for this. First, since the initial velocity is applied only to the 
particle and not to the interior fluid also, translational momentum MfU{0) and rotational 
momentum Jjfi(O) are quickly lost to the internal fluid which thereafter moves essentially 
as a rigid body. Second, the theoretical results apply to an incompressible fluid; thus there 
is the well known "added-mass" effect which accounts for the momentum carried off by 
sound waves in a real fluid. Our lattice-Boltzmann simulations are slightly compressible; 
thus a fraction of the initial momentum (MfU{0)/2) is carried off almost instantaneously, 
after which the velocity tracks the renormalized incompressible theory. Of course, there is 
no added mass for the rotational velocity. Thus the overall normalizations at t = are 
M/(M + 3M f /2)U(0) and 1/(1 + i/)O(0). 

VI. FLUCTUATIONS 

Suspensions of sub-micron sized particles undergo Brownian motion, due to the thermal 
fluctuations in the fluid. It has been shown (Ladd, 1993c), and further discussed in paper 
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I, that fluctuations can be incorporated into lattice-Boltzmann simulations by the straight- 
forward addition ol a random component to the fluid stress tensor. Thus the post-collision 
momentum flux in Eq 2.4 has in addition a fluctuating stress, <r', sampled from a Gaussian 
distribution and uncorrelated in space and time (Landau and Lifshitz, 1959), 

(cr^j(r, t)a' lS (r', t')^ = A6 rr >6 tt > {8 ai 8 rjS + 8 aS 8 rjl - (2/3)8 al3 8 lS ) ; (6.1) 

the variance A is proportional to the temperature of the fluctuating fluid. The exact re- 
lationship between the coefficient A and the effective temperature is discussed in paper I 
(Eq. 4.18). 

As a first test of the fluctuating lattice-Boltzmann simulations, the motion of an isolated 
sphere has been studied, at sufficiently short times that the periodic boundary conditions 
have no noticeable effect on the dynamics. The velocity correlation function of a spher- 
ical particle decays algebraically, with an asymptotic t~ 3 / 2 dependence which arises from 
hydrodynamic correlations in the fluid (Alder and Wainwright, 1970); here the fluctuating 
lattice-Boltzmann equation is tested to make sure that it can account for the hydrodynamic 
memory effects that lead to long-time tails. In Fig. 11 it can be seen that, within the sta- 
tistical error bars, the normalized velocity correlation functions are identical to the steady 
decay of the translational and rotational velocities of the sphere; thus our simulations satisfy 
the fluctuation-dissipation theorem. Once again there are no adjustable parameters in these 
comparisons. 

Next we consider the motion of suspensions of colloidal particles at very short times, prior 
to the onset of Brownian motion. Diffusing- wave spectroscopy (Zhu et al., 1992) has shown 
that if the mean-square displacement is normalized by the self- diffusion coefficient times 
the time (AR 2 (t)) /6D s (<f))t, and plotted vs. a reduced time t/r, then, for all solids volume 
fractions, there is a scaling time r(</>) which collapses the experimental data onto one master 
curve, indistinguishable from the isolated-sphere result. Moreover r seems related to the 
time it takes fluid vorticity to diffuse over a particle radius; values of r determined from the 
scaled mean-square displacements are in good agreement with independent estimates of the 
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vortex diffusion time ?y 00 (^)/ / oa 2 , based on the high-frequency suspension viscosity i]^. We 
have discussed this scaling in an earlier work (Ladd, f993c); here we just point out that the 
simulations reproduce the experimentally observed scaling over the whole time and density 
range. Simulation data for a number of Brownian particles (N = 128) under the same scaling 
is shown in Fig. 12. Results with different size systems (N = 16 and N = 1024) indicate 
that the periodic boundary conditions have a negligible effect on the 128-sphere results 
for times up to about 100r. The scaled data at various volume fractions collapse onto 
the dilute (single-particle) result, in excellent agreement with experiment both at very short 
times (Kao et al., 1993) and at somewhat longer times (Zhu et al., 1992). Moreover, the self- 
diffusion coefficient and viscosity that are required to scale the mean-square displacement are 
in quantitative agreement with independent simulations and experimental data (see Ladd, 
1990). A comparison is shown in Fig. 13. 

When the fluctuating lattice-Boltzmann equation is applied to a system of solid parti- 
cles, it is found that the translational and rotational velocities of the particles come into 
approximate thermal equilibrium with the fluctuations in the fluid. Because our model does 
not conserve the total energy, but instead preserves a balance, on average, between viscous 
and fluctuating forces, there is no equipartition principle to force an even division of energy 
between all the modes in the system. Thus we do not find an exact thermal equilibrium 
between the characteristic temperature of the fluid and solid motions (see Table VI). In ad- 
dition, the average kinetic energy of the solid particles is weakly dependent on the particle 
size and solids concentration; we find in general that the temperature characterizing trans- 
lation and rotation are similar, and that they are typically 10-20% less than the effective 
temperature of the fluid fluctuations. For studies of particle diffusion, temperature will be 
taken to be defined by the mean translational kinetic energy; for the shear viscosity, there 
is some ambiguity in defining the temperature which will be discussed later. 

The transport coefficients of a suspension of solid particles can be determined from 
the Green-Kubo relations (Hansen and McDonald, 1986). The simplest is the self- diffusion 
coefficient, which can be computed from the average mean-square displacement AR, or from 
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the velocity autocorrelation function; 

D s = lim 1 (AR(t) • AR(i)) = \ f°° (V(t) • U(0)) dt. (6.2) 

For comparison purposes, it is simpler to look at normalized correlation functions 

(U(f)-U(O)) 

/dW = (u(o).u(o))' (6 - 3) 

where macroscopically (U(0) • U(0)) = 3ksT/M; this expression is used to define the tem- 
perature of the suspension. A relaxation time td can be defined as the integral of the 
normalized velocity correlation function 

TD = / f D (t); (6.4) 

Jo 

thus the ratio D s /D = td/t 0} with r = M / (GTrrja). The integral in Eq. 6.4 and similar 
fluctuation expressions should be interpreted as a discrete quadrature; i.e. 

/ f D (t)dt = f D (0) + 2j2f D (2t), (6.5) 
Jo t=i 

where the factor of 2 arises because particle velocities are updated every two time steps. 

The ratio D s /D is shown in Fig. 14 for various volume fractions and particle sizes, using 
systems of 16 and 128 spheres. The agreement is quite good, but, compared with the steady 
non-equilibrium flows (Fig. 4), larger particles are required at high densities for the same 
accuracy. This implies that fluctuation-dissipation is not exactly satisfied at high solids 
concentrations (otherwise both measures of the diffusion coefficient would be the same); I 
suspect these discrepancies come from the effects of the shared nodes (section HID), which 
may well behave differently for fluctuating fluids than for steady flows. However, there is 
good agreement for the larger spheres, where the contribution from the shared nodes is 
negligible. 

The collective or mutual diffusion coefficient D is related to fluctuations in the total 
velocity of the solid phase 

N 

U C = £U,-. (6.6) 

8 = 1 
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The relationship between the fluctuations in U c and the diffusion coefficient D is discussed 
in the Appendix. Here we just quote the final result used for comparison with creeping-flow 
theory, 



where Xj is the mass fraction of the fluid. The results (Fig. 14) are quite comparable to 
the earlier results for collective mobility (Fig. 4) even at high concentration. This is further 
evidence that the relatively slow convergence of the self- diffusion coefficients in Fig. 14 arises 
from the shared nodes, which they play a minor role in collective diffusion because of the 
absence of lubrication. 

In section IV of paper I, it was shown that a discretized Green-Kubo relation could be 
derived relating the viscosity of the pure fluid to the equilibrium stress fluctuations; in our 
present notation, Eq. 4.13 of paper I can be written as 



where S-' is the fluid stress tensor. Equation 6.8 can be applied to solid-fluid suspensions 
by including both the fluid stress tensor (Eq. 4.12 of paper I) and the solid-fluid stress 
S s , summed over all the solid-particle surfaces. This latter contribution is computed in a 
similar way to the particle torques, summing the symmetric components of r&f (r;,) over all the 
boundary nodes. The fluid stress correlation function and the total stress (S* = S-^ + S s ) 
correlation function have been computed; by integrating these correlation functions the 
viscosity ratio ^(c^)/?? can be computed (Eq. 6.8). However, the estimated high-frequency 
viscosities for the suspension are then consistently too low, by up to 20% at high volume 
fraction. The reason is that energy is not uniformly distributed between the particles and 
the fluid, so that the effective temperature in Eq. 6.8 is different for the particle contributions 
and the fluid contributions. We have corrected for this by using the fluid temperature for 
the fluid-fluid contribution and the particle temperature (see Table VI) for the particle- 
fluid terms. A comparison of viscosities is shown in Fig. 15. Here we have averaged over 




(6.7) 
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all 5 components of the shear stress; in this case these angle- averaged viscosities are very 
insensitive to the number of particles so that results for 16 particles and 128 particles are 
indistinguishable. The agreement is again good, although the method of computation is 
not totally satisfactory. It may transpire that viscosities can only be computed reliably by 
some kind of external flow simulation as in section IIIC, and thus we would have to reserve 
the Green-Kubo method for quantities, like diffusion, which involve only particle-particle 
correlation functions. The question of the thermal equilibrium between solid and fluid needs 
to be researched in more detail; hopefully a method of including fluctuations will be found 
which succeeds in bringing about a true equipartition of energy between solid and fluid. 

VII. CONCLUSIONS 

The combination of molecular dynamics for the particulate phase and the fluctuating 
lattice-Boltzmann equation for the fluid phase has been shown to be a viable technique 
for quantitative simulations of hydrodynamically interacting particles. Even using small 
solid particles, with radii less than 5 lattice spacings, accurate results for hydrodynamic 
transport coefficients (permeability, sedimentation velocity, self- diffusion coefficient and vis- 
cosity) have been obtained over the whole range of packing fractions. The method is also 
very flexible; the particle size and shape, the electrostatic interactions, the flow geometry, 
the Peclet number and the Reynolds number, can all be varied independently. The technique 
is valid, without modification, at finite Reynolds number; lattice-Boltzmann simulations of 
flow past a column of cylinders are in quantitative agreement with finite-difference solutions 
at Reynolds numbers up to 100. The results for a dense array of spheres are in excellent 
agreement with finite-element calculations, but disagree with the truncated-mode analysis 
of Lahbabi and Chang. The simulations suggest that fluid inertia has a much smaller effect 
in periodic arrays than in random media. 

It has been shown that fluctuations can be incorporated into the lattice-Boltzmann 
equation in a very straightforward fashion, and then used to simulate the dynamics of 
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colloidal particles in suspension. The short-time transport coefficients of these suspensions 
can be accurately determined from autocorrelation functions of the fluctuations. We plan 
to examine long-time transport coefficients in the future, by tracking the motion of particles 
moving under the influence of Brownian forces. 
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TABLES 

TABLE I. Labeling of the 18 velocity model. Lor each label i, the velocity vector and speed 
(in lattice units) are shown. 
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TABLE II. Hydrodynamic radii of near- spherical objects. The friction coefficient ^ T = F^/Uy 
of a simple-cubic array of particles has been used to determine the effective hydrodynamic radius 
a in a suspending fluid with kinematic viscosity v = 1/6. The radius of the equal- volume sphere 
ay is also shown; the input radius <2q defines the boundary surface of the object. 



«0 a V 
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1.5 1.66 
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32 
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TABLE III. Hydrodynamic radii of near- 


-spherical objects as a function 


of fluid viscosity. The 


hydrodynamic radius of an object with input radius cio = 4.5 (ay 


= 4.53 


) is shown for various 


values of the kinematic viscosity v of the fluid. Numerical values 


for a were determined for a 


simple cubic lattice as in Table II. 








v a 
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1/24 4.67 
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1/12 4.60 
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TABLE IV. Steady flow past a column of cylinders up to Reynolds numbers R e = 100 
(R e = 2Uoa/i/). Simulations (LBE) for cylinders of radii a = 2.1, a = 5.1, and a = 10.1 are 
compared with accurate finite difference solutions (FD) of the time-independent Navier-Stokes 
equations (Fornberg, 1991) for the same channel width (or cylinder separation) W = 10a. In ad- 
dition to the drag coefficient C'd = -Fb/ pU^a, the wake length L\y, the minimum stream function 
$ c , and its location [X C ,Y C ], were determined from the contour plots (Fig. 6) of the recirculation 
zone behind the cylinder. At a Reynolds number R e = 100 the flow is eventually unsteady; the 
Strouhal number S = uj c a/Uo is also reported (uj c is vortex shedding frequency). 
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TABLE V. Steady flow past a simple-cubic array of spheres at maximum packing e = 1 — 7r/6. 



The volume- 


■averaged flow velocity Uy (Eq. 3.4) was 


computed for various pressure g 


;radients and 


sphere sizes. 


Results labeled by a 


drag force Fd = 
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creeping flow. 


The R e = 


friction coefficient (Zick and Homsy, 1982), F£>/(6in]aUv) = 42.1, can 
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The Reynolds number R e = 2Uva/ ' ev. 
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10.5 


33 


42.8 





33 


44.4 


10.5 


10 9 


50.5 


18.4 


20 9 


52.4 


35.4 


17 


45.9 


20.2 


17 


48.3 


38.4 


33 


46.6 


19.9 


33 


49.3 


37.6 


40 9 


55.2 


67.2 


70 9 


58.4 


111.2 


17 


50.9 


72.9 


17 


52.9 


122.9 


33 


52.2 


71.1 


33 


54.6 


118.9 


65 


52.5 


70.7 


65 


55.0 


118.1 


100 33 


56.4 


164.5 








65 


56.8 


163.3 


130 65 


58.4 


206.7 
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TABLE VI. Translational and rotational temperatures of suspended spheres, based on mean 
square translational and rotational velocities; the effective temperature of the fluctuating fluid is 



0.25. 



4> 


a 


N = 16 

T 


-*- T Ot 


N = 128 

T 




0.05 


1.54 


0.196 


0.164 


0.202 


0.166 




2.61 


0.215 


0.202 


0.218 


0.201 


0.25 


1.54 






0.183 


0.158 




2.61 


0.195 


0.196 


0.204 


0.196 




4.53 


0.202 


0.213 






0.45 


2.61 






0.183 


0.180 




4.53 


0.187 


0.206 


0.193 


0.203 




8.47 


0.207 


0.226 
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FIGURES 



Fig. 1. Geometry for quasi-periodic simulations. The figure illustrates a system that 
closely approximates a periodic one under an external flow. The lattice nodes are shown by 
solid circles and the boundary nodes by solid squares. The arrows indicate velocity directions 
Ci and c' at the boundary nodes (c.f. Fig. 2 of paper I). Periodic boundary conditions are 
applied across the planes indicated by solid lines; the configuration of solid particles within 
the quasi-periodic unit cells (bounded by dashed lines) are identical. Macroscopic flows can 
be set up by the planes of boundary nodes at either end the the system. With this geometry 
we can set up uniform flow perpendicular to the boundary walls or an approximately linear 
shear flow parallel to the boundary walls. The properties of the central cell are close to those 
of a truly periodic system; it is not necessary, apparently, to include more cells, although 
this can be done if required. 

Fig. 2. Translational and rotational friction coefficients of a simple-cubic lattice of 
spheres. The drag coefficients, normalized by the isolated sphere values, are plotted as a 
function of volume fraction for several different size objects. The solid lines are determined 
from accurate numerical solutions of the Stokes equations (Ladd, 1988; Ladd, 1990), as 
indicated in the text. 

Fig. 3. Hydrodynamic interactions between pairs of spheres. The perpendicular and 
parallel friction coefficients are plotted as a function of s = R/a — 2. The results at s = 
are for objects in closest possible proximity. The systems are periodic, with a two-sphere 
unit cell. The solid lines are again solutions of the Stokes equations in the same periodic 
geometry. 
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Fig. 4. Hydrodynamic transport coefficients for random arrays of spheres. Results from 
simulations of 16 spheres (with periodic boundary conditions) are compared with accurate 
numerical solutions of the Stokes equations (Ladd, 1990). The lattice-Boltzmann results 
(present work) are plotted as symbols; results from Ladd, 1990 are shown as solid lines. The 
statistical errors in both sets of calculations are smaller than the plotting symbols. 

Fig. 5. Velocity profiles in sheared suspensions at 5% and 45% solids concentration. The 
circles show the ensemble- averaged velocity profile; the solid circles indicate the central cell 
from which the viscosities shown in Fig. 4 are calculated. The velocities are scaled by the 
shear rate 70 of a sample of pure fluid under the same conditions. In dense suspensions, 
there is a rapid change in effective viscosity from the pure fluid regions (x < and x > 3L) 
to the suspension; thus the velocity profile is non-linear in this region. However the region 
surrounding the central cell has a linear velocity profile at all densities, although the local 
shear rate is less than 70. 

Fig. 6. Stream lines for steady incompressible flow past a column of circular cylinders 
at various Reynold's numbers. A small portion (1/8) of the total channel length, near the 
cylinder, is shown. 

Fig. 7. Pressure drop as a function of Reynolds number for steady flow past a simple- 
cubic array of spheres. The simulation results for different sized spheres are shown as solid 
symbols; the open symbols are the finite element results (Tompson, 1983). The solid line 
was fitted to the data of Lahbabi and Chang, 1985. 

Fig. 8. Frequency-dependent torque on a rotating sphere. The modulus and phase lag 
(in degrees) of the drag torque on a sphere undergoing small amplitude rotational oscillations 
are shown as a function of the reduced frequency coa 2 / 'v. Results are shown for two different 
size spheres at the same fluid viscosity v = 1/6; the solid lines are taken from expressions 
give in Landau and Lifshitz, 1959. 



36 



Fig. 9. Frequency-dependent torque on a rotating sphere as a function of the relative 
viscosity inside and outside the sphere. The modulus and phase lag (in degrees) of the drag 
torque on a sphere undergoing small amplitude rotational oscillations are shown as a function 
of the reduced frequency coa 2 ' jv exil for various viscosity ratios v r = Vmt/vext- Results are 
shown for viscosity ratios of 4, and 16; the hydrodynamic radii of the spheres are (from Table 
III) 4.67, and 4.79 respectively. The dashed lines include the full interior torque, whereas 
the solid line includes just the inertial contribution, corresponding to v r = oo. Note that 
the reduced frequency is based on the exterior fluid viscosity v exi . 

Fig. 10. Translational velocity U(t) and rotational velocity 0(t) of an isolated sphere. 
The time-dependent velocities of the sphere are shown as solid symbols; the solid lines are 
theoretical results, obtained by an inverse Laplace transform of the frequency dependent 
friction coefficients (Hauge and Martin-L6f, 1973) of a sphere of appropriate size (a = 2.61) 
and effective mass (p s / p = 11); the kinematic viscosity of the pure fluid v = 1/6. 

Fig. 11. Translational velocity correlation function (U(t)U(O)) and rotational velocity 
correlation function (O(t)O(O)) of an isolated sphere, normalized to their t = values. 
The time-dependent correlation functions are shown as solid symbols (with statistical error 
bars). The solid lines are theoretical results, obtained by an inverse Laplace transform of 
the frequency dependent friction coefficients (Hauge and Martin-L6f, 1973) of a sphere of 
appropriate size (a = 2.61) and mass (p s / p = 11); the kinematic viscosity of the pure fluid 
v = 1/6. 

Fig. 12. Scaled mean-square displacement (AR 2 (t)) /6D s t at short times, vs. reduced 
time t/r. Simulation results for 128 spheres (solid symbols) are shown at packing fractions <f> 
of 5%, 25% and 45%; the solid line is the isolated sphere result. The simulation parameters 
were the same as in Fig. 11 except that a sphere of radius 4.5 was used at the highest volume 
fraction. 
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Fig. 13. Scaled relaxation time t/t (t = a 2 jv) and self- diffusion coefficient D s /D 
vs. packing fraction <f>. The data points are determined from the scaling of the 128-sphere 
simulation data for various particle sizes. The uncertainty in fitting the data is about 5%. 
Independent results (Ladd, 1990) for rj /r] (solid lines) and D s /D (dashed lines) are shown 
for comparison. 

Fig. 14. Diffusion coefficients for random arrays of spheres. Self-diffusion coefficient D s 
and collective mobility // from simulations of 16 spheres (filled symbols) and 128 spheres 
(open symbols) are compared with numerical solutions of the Stokes equations (Ladd, 1990). 
Results from moment expansions of the Oseen equation are shown as solid lines (N = 16) 
and dashed lines (N = 128). The statistical errors in the fluctuating lattice-Boltzmann 
simulations are comparable to the size of the plotting symbols. 

Fig. 15. High-frequency shear viscosity for random arrays of spheres. Shear viscosities 
from simulations of 16 spheres (filled symbols) and 128 spheres (open symbols) are compared 
with numerical solutions of the Stokes equations (Ladd, 1990). The statistical errors in the 
fluctuating lattice-Boltzmann simulations are comparable to the size of the plotting symbols. 
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APPENDIX: COLLECTIVE DIFFUSION COEFFICIENT 



The Green-Kubo expression for the collective diffusion coefficient of a binary mixture can 
be derived from the fluctuations in concentration; here we outline the standard procedure 
given, for example, in Hansen and McDonald, f 986. The macroscopic equations of continuity 
for the two species, 

d t pi + V • ( Pt u t ) = 0, (Af ) 

can be cast into the form of an equation for conservation of mass (p = p\ + p 2 ) of the usual 
form, and an equation for the (mass weighted) concentration c = p\/p 

3 t c + u- Vc= -p-'V-y, (A2) 

j = ji = p\(ui — u) is the diffusive flux of species f, defined relative to the local stream 
velocity u = (p\Ui + p 2 u 2 )/p. The diffusive flux can also be written in a more symmetric 
form j = (f — c^iUi — cp 2 u 2 . 

A microscopic expression (in terms of molecular coordinates) for the spatial Fourier 
transform of the concentration can derived by linearizing the fluctuations in c about the 
mean densities ~p~i } 



r,.., P56pi(k,t)-pj8p 2 (k,t) X 2 E,- el e-*-*M 1 - X t E,- e2 e~^M 2 

oc(k,t) = == = = ; (A3) 

P P 



here Xi indicate average concentrations pj p, that are spatially invariant. Similarly, the 
diffusive flux can be expressed microscopically as 

j(M) = j e- tk ^(r } t)dr = X 2 J2e~ tk - R 'M 1 \J i - X 1 J2e~ tk - R 'M 2 \J i . (A4) 

Using the constitutive equation 

j = -pDVc (A5) 

to define the collective diffusion coefficient D, we find after some elementary manipulations 
that 
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/„" (J(f) • J(0)) <ft 

W (c(k)c(-k)) k=0 1 ] 

where J is the k = limit of Eq. A4 and the small k limit is taken in the denominator also. 
The equal-time concentration fluctuations can be expressed in terms of partial structure 
factors 

S ab (k) = . 1 y y e~ tk ^ ; (A7) 



p 2 (c(k)c(-k)) k=0 = X 1 X 2 M 1 M 2 In^O) - 2^nJj~ 2 S 12 (0) + N^O) 



(A8) 



The equal-time fluctuations in the diffusive flux are readily evaluated from the Maxwell- 
Boltzmann distribution of velocities, 

(J(0) • J(0)) = 3^TX 1 X 2 (iV 1 M 1 + N 2 M 2 ). (A9) 

We now specialize to the case of a particle suspension, for which the number of solid 
spheres (A^i) is much less than the number of solvent molecules (N 2 ); thus Eq. A8 simplifies 
to 

p 2 (c(k)c(-k)) k=0 = X 1 X 2 M 1 M 2 iV 2 5' 11 (0), (AfO) 

where Sii(k) is the structure factor for the solid particles. From Eqs. A6, A9, and AfO, 

k B T too (J(t) • J(0)) 
U ~ M 1 X 2 S 11 (0) Jo (J(0) • J(0)) M - 

In the simulations, fluctuations in J can be related to fluctuations in the velocity of the 
solid phase (by momentum conservation); i.e. J = Mi J2iei U 4 - = MiU c (Eq. 6.6). Using the 
relation 

/« = ^ (A12) 

k B l 

we arrive at Eq. 6.7 for the ratio p/ po. 



40 



The factor X 2 ~ X arises from the conservation of total momentum, which is implied in the 
definition of concentration fluctuations (Eq. A3). In a canonical ensemble, where the total 
momentum and energy are free to fluctuate, 

<U C (0) ■ U c (0)) = (A13) 

On the other hand, at fixed total momentum P, fluctuations in U c are reduced (Lebowitz 
et al., 1967), 

(U c (0) • U c (0)) p=o = (U c (0) • U c (0)) - (PP) : <9 P (U c ) p • d F (U c ) p 

= (U c (0) • U c (0)) - Mk B Tdp (U c ) p : d F (U c ) p , (A14) 

where Ai = N\M\ + N 2 M 2 is the total mass of solid and fluid. For fixed total momentum, 
the average particle velocity is ~P/ A4, and 

<9p(U c ) p = ^1. (A15) 

Combining Eqs. A16 and A15, the fluctuations at constant total momentum are found to 
be 

(U,(0).U,(0)) p=0 = 5^Zl. (A16) 

Thus the ensemble dependence of the zero-time fluctuations accounts for the factor X 2 X \ in 
an ensemble where the total momentum is free to fluctuate this factor is absent. 
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